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Abstract. The magnetization switching dynamics in the kinetic Ising model is 
projected onto a one-dimensional absorbing Markov chain. The resulting pro- 
jected dynamics reproduces the direct simulation results with great accuracy. 
A scheme is proposed to utilize simulation data for small systems to obtain the 
metastable lifetime for large systems and/or for very weak magnetic fields, for 
which direct simulation is not feasible. 

In simulations of metastable decay one faces the problem of measuring the life- 
time of the metastable phase, which is by definition very long. For example, 
in Monte Carlo modeling of magnetization switching in ferromagnets the phys- 
ically relevant simulation time scales are on the order of 10 12 — 10 15 Monte 
Carlo Steps per Spin (MCSS). Even with sophisticated algorithms [1,2] such 
simulations have not yet been feasible, and one has to resort to extrapolation 
of the results into the physical time-scale regime. 

In this article, a scheme is presented to map the system under study onto a 
simpler one, which is faster to simulate but still gives accurate results for most 
important physical quantities. This idea is not new. For example, recently 
Lee et al. introduced a macroscopic mean-field dynamics [3] which semiquan- 
titatively describes magnetization switching in Ising systems. The aim of the 
present work is to construct a scheme which is very similar in spirit, but is 
intended to create a practical computational tool applicable to the simulation 
of ferromagnets in the physical time regime. 

To simplify our notation and reasoning, we explain the proposed Projected 
Dynamics (PD) as applied to the isotropic kinetic Ising model on a square or 
cubic lattice. The Hamiltonian has the standard form, 



with the ferromagnetic nearest-neighbor spin-spin interaction J > and an 
external magnetic field H. In what follows, we denote by V the total number of 
spins in the system. To study the magnetization reversal, we initialize all spins 
in the state +1, fix the temperature T well below its critical value T c , and apply 




(1) 



a negative magnetic field. Then we apply Metropolis or Glauber dynamics with 
updates at randomly chosen sites to measure the time the system needs to reach 
a configuration with a given stopping magnetization. Repeating this procedure 
many times, we obtain the mean lifetime, t, of the metastable state and its 
standard deviation, At. 

To speed up the simulations as much as possible, we use a rejection-free 
algorithm [1,2] which uses the notion of spin classes. By a spin class we mean 
the state of the spin itself and its neighbors. There are ten classes for the 
given model on a square lattice. Classes i = 1, . . . , 5 correspond to spins in the 
state +1 which have exactly i — 1 neighbors in the state +1. Similarly classes 
i = 6, . . . , 10 are assigned to those —1 spins which have i — 6 neighbors in the 
state +1. All spins in class i have the same flipping probability, pi. 

Rejection-free algorithms keep track of the number Cj of spins in each class, 
and it is not computationally expensive to measure the growth and shrinkage 
rates of the stable phase, which are defined as 
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9(n) = ^2{ci) n Pi , s(n) = ^2(ci) n pi , (2) 
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respectively. The angular brackets mean the average taken over the configura- 
tions generated during the Monte Carlo lifetime experiment, conditionally on 
the number n of overturned spins. Thus, g(n)/V is the probability that a +1 
spin will be flipped in the next Monte Carlo step, and s(n)/V is the probabil- 
ity to flip one of the —1 spins, both conditionally on n. Flipping a +1 spin 
increases the volume fraction of the stable phase, hence the names g(n) and 
s(n). 

The main idea of the proposed method is to make use of the observed 
growth and shrinkage rates to map the switching dynamics onto a one-di- 
mensional absorbing Markov chain. We assign to all configurations with n 
overturned spins a single state n in the chain. The one-dimensional dynamics 
is given by the above probabilities: From the state n we have the probability 
g(n)/V of jumping to the state n + 1, the probability s(n)/V of jumping to 
n — 1, and the probability 1— g(n)/V — s(n)/V of remaining in the current state. 
This random walk starts at n = and terminates when it reaches n = N + 1 
where M = V — 2N — 2 corresponds to the stopping magnetization. Using 
standard methods from the theory of absorbing Markov chains [2,3], we obtain 
the mean lifetime r and the total average time h(n) (in MCSS) spent by the 
random walker in the state n in terms of the growth and shrinkage rates: 

E\_ , , . 1 + s(n + l)h(n + 1) , 1 
, =n • h{n} = P(n) ■ w-m- (3> 



Higher moments of the lifetime distribution can be obtained in a similar way. 
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Fig. 1 Illustration of the Projected Dynamics [Eqs. (2,3)] and of the "size 
extrapolation" scheme [Eq. (4)]. The lines in both panels are Projected Dy- 
namics extrapolations based on (cj)„ sampled during lifetime measurements on 
a 20 x 20 lattice. Symbols are direct simulation results, a) The lifetime and 
its relative standard deviation vs. the system volume. The agreement is near- 
perfect in the single-droplet regime (where Ar/r « 1), and deviations are only 
observed in the crossover region to the multidroplet regime (where At/t — > 0). 
b) Extrapolation from the 20 x 20 lattice reproduces the measured growth and 
shrinkage rates, g(n) and s(n), of the 160 x 160 lattice very well. The crossings 
of the two curves correspond to the metastable magnetization (left) and the 
critical fluctuation (right). 

As shown in Fig. 1, the proposed scheme reliably reproduces the direct 
Monte Carlo results. The lifetimes obtained from Eq.(3) are usually only 
slightly different from their counterparts from the direct measurements. The 
projected dynamics also gives the standard deviation of the lifetime distribu- 
tion, which agrees with the measured values within the error bars. 

Suppose we have measured g(n) and s(n) in a system of volume V. How are 
they related to their counterparts in a system of volume 2V? Since the relevant 
configurations typically contain many small droplets of the stable phase, it is 
a reasonable approximation to view the larger system as consisting of two 
"independent" copies of the smaller system. Then, the growth rate of the large 
system can be approximated as a weighted sum of the sub-system contributions, 



g(2V,n) 



£?=o HV, n - i)h(V, i){g(V, n-i)+ g(V, i)} 
EtoHV,n-i)h(V,i) 



(4) 



Here, we have explicitly shown the dependence on the volume V. An analogous 
formula is proposed for the shrinkage rate s(n). Provided the smallest system 
is in the so-called single-droplet regime, in which the nuclcation is triggered 
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Fig. 2 Populations of two spin classes vs. the number of overturned spins 
n. Solid lines are (ci) njG q measured in the equilibrium fixed-7z ensemble; circles 
and dotted lines show (ci) n measured during the lifetime measurement for two 
different fields. Note that as the field decreases the region in which (c,)„ — ► 
(c%)n,eq increases. 

by a single critical fluctuation smaller than the system size (for the theoretical 
description of different regimes of the magnetization switching, see Refs. [4- 
8]), one can repeat the extrapolation V — > 2V several times without finding a 
big discrepancy between growth and shrinkage rates calculated from the small- 
lattice data and those directly measured on large lattices (see Fig. lb). 

Now we return to our main goal, namely predicting lifetimes in very weak 
fields. What we need is to calculate g(n) and s(n) as functions of the magnetic 
field. The naive choice would be to replace the class spin-flip probabilities 
Pi by their values calculated for the desired field, and keep the values (cj)„ 
entering Eq. (2) unchanged. Such a scheme works quite well when extrapo- 
lating to stronger fields, but the lifetimes are systematically overestimated for 
fields weaker than the one at which the (c,)„ were sampled. The reason lies 
in the noncquilibrium character of the configurations with n larger than the 
critical droplet volume (see Fig. 2). As it tries to escape from the metastable 
free-energy minimum, when n is small the system passes through configura- 
tions which are very close to "equilibrium" . Farther from the free-energy min- 
imum, the configurations which appear in the system are increasingly different 
from "equilibrium" configurations. Although these configurations do not con- 
tribute significantly to the lifetime, if we use them to estimate the lifetime in 
a weaker field their contribution becomes more important because the system 
then spends more time in them. However, the actual configurations are closer 
to "equilibrium," because the free-energy minimum is deeper in a weaker field. 




30 



d=2, T=1.3J=0.57T C 
Metropolis dynamics 



d=3, T=2.0J=0.44T C 
Glauber dynamics 



30 



20 



10 











0.0 



2.0 



4.0 



.0 6.0 
1/H [J 1 ] 



8.0 



10.00.0 



0.5 



1.0 

1/H [J ] 



1.5 



Fig. 3 The metastable lifetime of a kinetic Ising model as a function of 
the magnetic field. Points are conventional Monte Carlo measurements, and 
lines are Projected Dynamics calculations based on the "equilibrium" (ci)„. oq 
data sampled on the smaller systems. Predictions of the Projected Dynamics 
improve with decreasing field. 

Thus, we effectively replace the actual weak-field configurations by ones which 
appear slightly "overheated." This enhances the shrinkage probabilities more 
than the growth probabilities, thus leading to an overestimate of the metastable 
lifetime. Understanding the cause of this problem also offers a remedy: we 
can do better by using the equilibrium configurations. Then we expect to see 
underestimation of relatively short lifetimes (by reversing the above argument). 
On the other hand, such an approximation will improve with decreasing field, 
which is exactly what we need to extrapolate towards zero field. Thus, we 
define the equilibrium growth and shrinkage rates as follows 



The only difference from Eq. (2) is in the sampling. Here, (ci)„ jeq are sampled 
in an equilibrium ensemble with a fixed number of overturned spins n. Thus, 
instead of sampling the configurations during the lifetime measurement, we 
perform a static measurement for each value of n needed (up to the stopping 
magnetization). The (ci) n ,eq depend on temperature, but the only dependence 
of g cq and s oq on the magnetic field comes from the spin-flip probabilities pi. 
In that way, a single measurement provides sufficient data to obtain a good 
approximation for the lifetime as a function of the external field. 

Figure 3 shows the lifetime calculated from the projected dynamics for 
two pairs of lattice sizes in two and three dimensions. Comparison with the 
direct simulation results (points) corroborate our expectation that the lifetime 
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is underestimated in the crossover region between the single-droplet and mul- 
tidroplct regimes, and that the approximation provides progressively better 
results as the external field strength decreases. Measurements on small lat- 
tices, for which we can obtain reasonable statistics in zero field, show that the 
Projected Dynamics based on (ci)„ jCq reproduces the "lifetime" in zero field. 

In conclusion, our Projected Dynamics maps the complex Monte Carlo 
dynamics onto a much simpler one-dimensional absorbing Markov chain. The 
Projected Dynamics is computationally much easier to study than the full un- 
derlying Monte Carlo dynamics, while it provides reliable results for metastable 
lifetimes and their standard deviations. Data from small systems can be uti- 
lized to predict lifetimes for large systems. When based on "equilibrium" class 
populations (ci)„ jCq , Projected Dynamics yields results as functions of the field. 
Such estimates are most reliable in the experimentally relevant weak-field re- 
gion, which it is currently not feasible to investigate with direct simulations. 
The method also offers a deeper insight into the mechanism of metastable de- 
cay. For example, the histogram h(n) provides a simple method to measure 
the metastable magnetization and susceptibility, and the rates g(n) and s(n) 
contain information about the size of the critical fluctuations. 

For other types of models, the applicability of the projected dynamics 
will depend on the existence of a "single channel" for the decay. The pro- 
jected parameter (n in the present case) need not necessarily be related to the 
magnetization, but it should parameterize the optimal path for escape from 
metastability in such a way that the configurations at a fixed value of the 
parameter will exhibit similar growth and shrinkage rates. 
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